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Abstract 

This  paper  deals  with  the  derivation  of  equations  suitable  for  the 
computation  of  elastic  curves  on  the  sphere.  To  this  end  equations 
for  the  main  invariants  of  spherical  elastic  curves  are  given.  A  new 
method  for  solving  geometrically  constraint  differential  equations  is 
used  to  compute  the  curves  for  given  initial  values.  A  classification  of 
the  fundamental  forms  of  the  curves  is  presented. 


1  Introduction 

Two  current  trends  in  Geometric  Modeling  are  concerned  with 

•  the  development  of  spline  techniques  on  surfaces  (see  [14],  [12],  [10]) 

•  the  use  of  nonlinear  spline  curves  of  minimal  elastic  energj'  for  the 
modeling  of  smooth  shapes  (see  [2],  [3],  [4]). 

This  paper  blends  both  topics  in  considering  elastic  curves  on  the  sphere. 

In  euclidean  space  an  elastic  curve  can  be  viewed  as  an  arc  length  para¬ 
metrized  "cubic”  spline  in  tension,  i.e.  an  elastic  curve  is  a  critical  point  of 
the  functional 

^(y)=f  <  y",y"  > <  y\y  >  ds  (1) 

Jo 
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in  the  space  F  of  smooth  maps 


y;[0,/]-^R",  |y'l  =  l 

y(0)  =  Po,  y(i)  =  Pi,  y  (0)  =  Vo,  y'(l)  =  Vi 

where  Pq,  Pi  €  R",  Vq  G  Tp^R",  Vi  G  TpjR”,  a  G  R  are  fixed  and  I  is 
variable.  Expressing  y"  in  the  Frenet  frame  yields  the  functional  $  in  the 
form 

=  I  {ti{s))^  +  a  ds  (2) 

Jo 

which  explains  the  interest  in  elastic  curves  as  those  curves  that  minimize 
bending. 

The  notion  of  cubic  splines  can  be  generalized  to  curves  on  a  Riemannian 
manifold  M  by  replacing  the  usual  derivative  of  the  tangent  vector  field 
y'  by  the  covariant  derivative  compatible  with  the  metric  of  M  (see  [13]). 
Generalizing  the  functional  (1)  in  this  way,  one  obtains  the  concept  of  elastic 
curves  on  arbitrary  manifolds.  In  the  case  of  surfaces  embedded  in  R^  the 
algebraic  value  of  the  covariant  derivative  of  the  tangent  vector  field  y'  of  a 
surface  curve  y  is  called  the  geodesic  curvature  iCg  of  y  (see  [5]  and  section 
2).  Therefore  we  may  define  an  elastic  curve  on  a  surface  5  :  /I  C  R^  — >  R^ 
as  an  extremal  point  of  the  functional 

♦(»)=  +  (3) 

Jo 

in  the  subset  P  of  P  formed  by  curves  lying  on  S. 

In  section  3  a  set  of  differential  equations  for  eleistic  curves  on  the  sphere  is 
derived.  This  set  includes  a  differential  equation  for  the  geodesic  curvature  of 
spherical  elaistica.  Since  the  normal  curvature  of  a  spherical  curve  is  constant, 
the  differential  equation  for  the  geodesic  curvature  suffices  to  compute  the 
ordinary  curvature  of  a  spherical  elastica.  Furthermore,  a  formula  is  given 
that  expresses  the  squared  torsion  of  a  spherical  elfistica  as  a  rational  function 
of  its  curvature. 

In  section  4  we  describe  the  numerical  algorithm  used  to  integrate  the 
set  of  differential  equations  derived  in  section  3.  The  equations  have  a  very 
particular  structure  defined  by  a  number  of  constants  of  motion,  and  in 
particular  they  constrain  the  elastic  curves  to  lie  on  the  sphere.  We  employ 
an  algorithm  introduced  by  Crouch  and  Grossman  [8]  which  preserves  the 
constraints  exactly. 

Since  Euler’s  fundamental  work  on  plane  eleistic  curves  it  is  known  that 
these  curves  can  be  classified  according  to  their  shape.  The  tools  developed 
in  this  paper  enable  us  to  present  the  fundamental  forms  of  spherical  elcistica 
in  the  Icist  section. 
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2  Geometric  Preliminaries 


Let  5  :  C  R2  denote  a  regular  parametric  surface  and  N  the  unit 

normal  vector  field  of  S.  A  curve  i  :  /  C  R  — ♦  R^  is  a  curve  on  the  surface 
S  if  and  only  if  x  =  5  o  c  where  c  :  I  —*  (J  is  a.  plane  curve  in  U.  The  unit 
normal  of  S  along  a  surface  curve  x  will  be  denoted  hy  n  :=  N  o  c. 

The  Darboux  frame  61,62?  ^>3  along  x  is  the  orthogonal  frame  defined  by 

x'(t) 

=  b,{t)  =  n{t)  X  b,{t),  63(<)  =  n(0. 

F  (i)l 

The  equations  that  express  the  derivatives  61,63,63  in  the  Darboux  basis 
61, 62,  63  are  given  by; 


61  =  UJKgb^  +  OJKnbs,  (4) 

62  =  —UJKgbl  +  UJTgb3,  (5) 

63  =  -UJKgbl -(jJTgbi  (6) 

with  w(f)  =  |3:'(<)1. 

The  functions  Kg,  and  Tg  are  called  geodesic  curvature,  normal  cur¬ 
vature  and  geodesic  torsion.  The  geodesic  torsion  and  the  absolute  value  of 
geodesic  and  normal  curvature  are  invariant  under  reparametrization  of  the 
surface. 

The  geodesic  curvature  of  a  surface  curve  i  at  a  point  x[t)  is  the  ordinary 
curvature  of  the  plane  curve  generated  by  orthogonal  projection  of  x  onto 
the  tangent  plane  of  5  at  i(t).  It  can  be  computed  using  the  formula 


Kg  = 


X  ,  X  ,  nj 


(7) 


A  surface  curve  with  identically  vanishing  geodesic  curvature  is  called  a 
geodesic  of  the  surface. 

The  absolute  value  of  the  normal  curvature  of  x  at  a  point  x{t)  is  the 
curvature  of  the  intersection  of  5  with  the  plane  through  x{t)  spanned  by 
the  vectors  x'{t)  and  n(t).  While  the  geodesic  curvature  is  the  curvature  of  a 
surface  curve  from  a  viewpoint  in  the  surface,  normal  curvature  measures  the 
curvature  of  the  curve  that  is  due  to  the  curvature  of  the  underlying  surface. 
If  K  denotes  the  ordinary  curvature  of  the  space  curve  x  the  identity 


(8) 


holds. 
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The  geodesic  torsion  of  a  surface  curve  x  at  a  point  x(t)  is  the  torsion 
of  the  geodesic  that  meets  x  at  x{t)  with  common  tangent  direction.  A 
curvature  line  of  x,  i.e.  a  curve  with  a  tangent  vector  that  points  into  one  of 
the  principal  directions  of  the  surface,  is  characterized  by  vanishing  geodesic 
torsion. 


3  The  differential  equations  of  spherical  elas- 
tica 

Let  5  be  a  parametrization  of  a  patch  on  a  sphere  S2  C  of  radius  r  and 
center  0  such  that 

S  =  rN 

and  X  be  an  arc  length  parametrized  (i.e.  ]x'|  =  1)  curve  on  S.  Then 

X  =  rn  = 

In  this  situation  (6)  implies  that 

Tg  =  0.  (9) 

r 

Since  the  absolute  value  of  and  Tg  are  invariant  under  reparametrization, 
these  equations  imply  that  any  spherical  curve  is  a  curvature  line  with  con¬ 
stant  normal  curvature. 

We  consider  the  variational  problem  of  minimizing 

f{Kg{s)r  +  ads 

Jo 

in  the  space  F  of  C°°  smooth  maps 

y  :  [0,/]  -♦  S2,  |t/'|  =  1 

y(0)  =  Fo,  y{i)  =  Pu  y'{0)  =  Vo,  y{\)  =  V 

where  Pq,  Pi  £  S2,  Vo  €  Tp,S2,  V’l  €  Tp^S2,  <7  €  R  are  fixed  and  I  is 
variable. 

From  (4)  one  obtains  the  relation 

=  + 

SO  that  we  wish  to  minimize  the  functional 

[‘ \b\{3)\^  +  S  ds 

Jo 
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where 


under  the  constraints 


S  —  a  — 


l&il^  =  l,  h  =  |xj^  =  r^. 

Hence,  we  can  apply  the  Euler- Lagrange  equations  to  the  functional 
F  =  |6;|2  +  5  +  A(|6i|2  -  1)  -h  -  r^)  +  2  <  A, x'-b,> 
to  obtain  the  differential  equations  which  govern  the  extremals: 


Hx  -  A'  =  0 

(10) 

A61  -  6"  =  A. 

(11) 

Combining  these  equations  yields 

X'bi  +  Xb[  —  b'"  =  fix. 

(12) 

The  derivatives  of  bi  expressed  in  the  Darboux  basis  are  given  by: 

K 

=  /C362  -  -63 

r 

(13) 

K 

(14) 

6r 

=  (  — +  (k"  —  —  '^'^g)^2 

(15) 

Substituting  these  derivatives  and  x  =  rb^  into  (12)  and  rearranging  gives 


(A'  -f  3KgK'^)bi 

-  (i  +  +  4))6,  =  0. 

r  r  r 

Finally,  the  linear  independence  of  the  vectors  6j,  62  and  63  implies  that 

a  =  -5«j  +  c  (16) 

and 

<  +  =  »■  ('■) 
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In  order  to  determine  the  constant  C  in  terms  of  the  tension  parameter 
<7,  we  consider  the  boundary  condition 

m  -  |^(0x'(/)  -  =  0 

for  the  extremal  x.  This  condition  is  implied  by  the  fact  that  the  total  length 
I  of  the  curve  is  variable  in  the  variation  (see,  e.g.,  [1]).  Thus, 


-  kIU)  -±  +  S-2<  A(/),x'(/)  >=  0.  (18) 

Substituting  A  according  to  (11)  into  the  scalar  product  <  A.i'  >  yields 

<A(/),x'(/)>  =  <xb,{i),b,{i)>-<b';{i),h{i)> 

=  A  +  («:^(/))^  +  ^ 

=  -|w'))’  +  ^  +  c- 

Substituting  this  expression  for  the  scalar  product  into  (18)  we  obtain 

These  results  are  summarized  in  the  following  theorem. 

Theorem  1  An  elastic  curve  x  under  tension  a  on  the  sphere  of  radius  r 
satifies  the  differential  equations 

^  I'l  \  /  0  l/r  0  \  /  xi  \ 

^2  =  -1/''  0  /Cj  I  X2  \  (20) 

\  x's  J  [  0  -Kg  0  I  \  X3  ) 

where  xi  =  x,  X2  =  rx' ,  x^  =  x  x  x'  and  where  the  geodesic  curvature  Kg  of 
X  is  a  solution  of 

'^3  +  +  (rr  “  (2^) 


The  curvature  of  a  spherical  elastic  curve  can  be  obtained  from  (8)  and  the 
fact  that  the  normal  curvature  of  spherical  curves  is  constant.  The  squared 
torsion  of  a  spherical  elastica  can  be  expressed  as  a  rational  function  of  its 
curvature. 

Theorem  2  The  curvature  k  and  the  torsion  t  of  a  spherical  elastic  curve 
obey  the  relation: 


h^K*  =  -  i(l  -  a)K^  +  ^(  J  + 
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Proof;  For  the  invariants  k  and  r  of  an  arc  length  parametrized  curve  x  in 
the  relation 

TK^  =  (22) 

holds  where  [a,b,c]  denotes  the  determinant  of  three  vectors  in  R\ 

For  a  spherical  curve  b\  and  b'(  can  be  expressed  in  the  Darboux  basis  as 
follows: 


K  = 
K  = 


Kgb2  - 


r 


+  +  k' ^2- 


Substituting  the  derivatives  into  (22)  and  squaring  yields 


Since  the  differential  equation  (21)  can  be  integrated  to 


(4- 


r2  2'"^’ 


one  obtains  the  claimed  equation  using  (8). 


4  Tracking  elastic  curves  on  the  sphere 

The  problem  w'e  consider  here  is  that  of  numerically  integrating  equations 
(20)  and  (21)  in  Theorem  1.  One  can  of  course  simply  integrate  the  equa¬ 
tions  using  a  standard  numerical  package,  such  as  an  IMSL  Runge  Kutta 
routine.  However,  the  system  of  equations  possesses  a  very  special  structure. 
As  pointed  out  in  [3],  equation  (21)  may  be  integrated  directly  in  terms  of 
Jacobi’s  elliptic  functions.  We  give  more  details  of  this  process  in  the  next 
section  where  we  classify  the  various  extremals.  As  for  equation  (20)  we 
note  that  the  components  of  the  state  vector  [x^, satisfy  algebraic 
constraints  consistent  with  the  fact  that  the  matrix  [zi,Z2,Z3]  is  simply  a 
multiple  r,  of  a  rotation  matrix.  When  a  standard  integration  package  is 
applied  to  the  set  of  differential  equations  (20)  and  (21),  these  constraints 
are  not  preserved  exactly,  and  in  particular  the  norm  of  the  vector  Zj  will 
not  remain  at  the  constant  value  r.  This  is  a  particularly  important  fact 
when  we  wish  to  integrate  the  equations  over  a  large  number  of  time  steps 
and  visualize  the  resulting  curves. 

We  have  therefore  made  use  of  a  new  class  of  integration  algorithms  de¬ 
veloped  by  Crouch  and  Grossman  [8],  and  Crouch,  Yan  and  Grossman  [7] 
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which  do  indeed  preserve  such  structures.  The  algorithms  are  therefore  called 
geometrically  exact  [6].  We  briefly  indicate  the  important  cispects  of  these 
geometrically  stable  algorithms  which  pertain  to  the  equations  (20)  and  (21). 
.Suppose  that  we  wish  to  numerically  integrate  an  ordinary  differential  equa¬ 
tion  on  R”  given  by  the  equations; 

i(<)  =  F(/,x(t)),  xGR\  t{0)  =  xo  (23) 

where  ^ 

F(t,x)  =  ^a-'(f,x)/lj(i)  (24) 

n  >  l\\  Aj  are  vector  fields  on  R"  and  are  functions  on  R  x  R".  Suppose 
that  in  addition  we  are  given  a  set  of  functions  on  R”  whose  numerical  values 
are  constant  along  solutions  of  the  equation  (23)  and  that  the  level  sets  of 
these  functions  are  manifolds.  Denote  the  level  set  through  Xo  by  M .  It  is 
convenient  to  assume  the  slightly  stronger  assumption  that  the  vector  fields 
Aj  are  everywhere  tangent  to  M.  We  also  assume  that  there  is  an  oracle  that 
can  integrate  any  vector  field  of  the  following  form,  to  any  desired  accuracy: 

Z(x)  =  f:aM,(x).  (25) 

j=i 

Here  q-'  are  real  numbers.  We  define  vector  fields  by  setting: 

F'(x)  =  f:aJ(p)/l,(x) 
r=i 

and  note  that  is  simply  the  vector  field  F  with  coefficients  “frozen”  at 
the  point  p.  If  we  denote  the  flow  of  any  vector  field  Z  by  (f,x)  —*  Bz[t^x)^ 
R  X  R”  R",  then  since  the  vector  fields  Aj  are  everywhere  tangent  to  A/, 
it  follows  that  x  (r  M  implies  that  9pp{t,x)  €  M  for  all  p  and  for  all  t  for 
which  the  flow  is  defined. 

We  now  introduce  the  (explicit)  geometrically  exact  Runge  Kutta  algo¬ 
rithms  as  described  in  [8].  Let  Xk  =  x{tk)  be  a  point  of  the  integral  curve 
X  of  (23).  Then,  define  vector  fields  on  R"  by  freezing  coefficients  of  F  at 
various  points  cis  follows: 

N 

j=i 

N 

+  hc2uOFi{hc2uXk))Aj{x) 

j=i 
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A' 

P'six)  =  +  C32),&F2(hc32,0FA^(^31iP)))Aj{x), 

J  =  1 

etc.  Second,  we  define  the  numerical  integration  algorithm  via  an  update 
rule; 

^k+\  =  0Fr{hCr-,0f^_^{hcr-i,-  ■  ■  ,eF^{hcuXk)))  (26) 

wlierc  h  is  the  “step  length”  and  and  Cjj  are  constants  to  be  determined. 
These  constants  are  determined  from  the  “consistencj'  equations,”  obtained 
by  making  the  Taylor  expansions  in  h,  about  /i  —  0,  of  both  sides  of  (26), 
using  on  the  left  hand  side  the  expression  xjt+i  =  6p{h,xi;).  If  the  coefficients 
of  h'  agree  up  to  i  =  q,  then  q  is  said  to  be  order  of  the  resulting  algorithm. 
Note  that  in  general  we  have  r  >  q,  while  for  classical  Runge  Kutta  schemes 
we  can  always  take  r  =  q.  Note  that  the  update  rule  defined  by  equation 
(26)  has  the  property  that  if  x*,.  lies  in  A/,  then  so  does  x^+i  since  each  flow  is 
defined  by  a  vector  field  F  with  frozen  coefficients.  In  the  special  case  where 
71  =  A’,  M  =  R"  and  A,  =  e,  is  the  standard  rth  basis  vector  in  R",  the 
algorithm  reduces  to  the  form  of  a  classical  explicit  Runge-Kutta  algorithm. 

In  the  paper  [8]  the  consistency  equations  are  derived  for  the  geometri¬ 
cally  exact  Runge  Kutta  algorithms  via  a  careful  geometric  analysis  of  the 
equations  (26).  The  results  show  that  a  geometrically  exact  third  order  ex¬ 
plicit  Runge  Kutta  algorithm  can  be  obtained  for  r  =  3  and  is  determined 
by  five  independent  consistency  equations,  in  the  six  constants  defining  the 
algorithm.  The  equations  have  multiple  solutions,  all  of  W'hich  are  solutions 
to  the  equations  which  determine  the  classical  explicit  Runge  Kutta  algo¬ 
rithms.  Thus,  all  solutions  define  classical  algorithms,  but  the  solutions  are 
not  ones  traditionally  found  in  the  literature.  We  have  used  the  following 
solution  of  all  five  equations: 

Cl  =  1,  C2  =  —2/3,  C3  =  2/3, 

C21  =  —1/24,  C31  =  161/24,  C32  =  —6. 

In  the  special  case  of  the  set  of  equations  (20)  and  (21),  defining  the 
elastica  on  a  sphere,  we  use  a  hybrid  algorithm  in  wdiich  equations  (20)  are 
integrated  using  the  third  order  geometrically  exact  Runge  Kutta  algorithm 
described  above  by  freezing  the  coefficients  Kg.  These  coefficients  are  then 
updated  by  integrating  equation  (21)  using  elliptic  functions. 

In  equation  (20),  M  is  the  three  dimensional  submanifold  of  deter¬ 
mined  by  the  equations 


<  xi,x,  >= 


<  X,,X2  >=  0, 


<  X2,X3  >=  0, 


<  X2,X2  >=  1, 
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<  Xi,X3  >=  0, 

<  X3,l3  >=  1. 


The  solutions  of  (20)  for  any  initial  condition  x(0)  =  Xq  G  A/  lie  completely 
in  M.  This  follows  from  the  fact  that  the  functions 


/l=<X,,Xi>,  /2  =<  >,  /3=<Xl,X3>. 

fs  >1  /e  =<  2r3,X3  > 

satisfy  the  system  of  differential  equations 

/;  =  v-2 

/2  =  /4  +  K3/3-{l/r^)/l 

fs  =  /s  “  s/2 
/'  =  2Kj5-(2/r^)/2 
/s  =  s/e-(lA'")/3-s/4 
/6  =  -2s/5 

which  has  the  unique  solution 

/,  =  r^  h  =  0,  /3  =  0,  /4  =  1,  /s  =  0,  /e  =  1. 

Since  this  argument  does  not  specify  the  function  Kg,  it  follows  that  the  flows 
of  the  vector  fields  F\  —  F3  with  frozen  coefficients  are  mappings  into  M  as 
needed  in  formula  (26). 

Furthermore,  the  vector  fields  Aj  in  equation  (20)  are  linear,  so  that  F  is 
given  by  an  expression  of  the  form 


F{x)^'£b’{t,x)BgX  (27) 

>=i 

where  Bj  are  matrices  and  are  functions.  Freezing  the  functions  b’  to  values 
/?■’  yields  a  system  of  linear  differential  equations  with  constant  coefficients: 

i=(f;^^5,)x.  (28) 

^3  =  1  ' 

Thus,  the  flow  of  the  vector  field  Fi  (i= 1,2,3)  is  given  by 

=  exp(fC,)9 

and  (26)  takes  the  special  form 


2Ja:+i  =  exp(c3/iC3)  •  exp(c2/iC'2)  ’  exp(ci/iCi)  ■  Xfc. 
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Since  6’  =  Ijr  is  constant  and  6^  =  Hg  depends  only  on  the  matrices  Ci 
are  given  by 

Cl  =  C{tk),  C2  —  C{tk  +  hc2l),  C3  =  C(tji:  +  h(C3i  +  C32)) 

where  C(t)  :=  5Ij=i  b^{t)Bj.  C  can  be  considered  as  a  3  x  3  skew  symmetric 
matrix  with  matrix  components  that  are  themselves  3x3  matrices.  Therefore 
the  standard  formula  for  the  exponential  of  a  3  x  3  skew  symmetric  matrix 
can  be  used  to  compute  the  flow  of  each  vector  field: 

exp(<05(c))  =  /  +  sin{<c>).S'(c)  +  (1  —  cos{t(p))S{c)^ 

where  5(c)  is  the  skew  symmetric  matrix  satisfying  S(a)b  =  bxa.  and  |c|  =  1. 

Regarding  the  performance  of  the  geometrically  stable  integration  method 
it  is  clear  that  this  method  is  computationally  more  expensive  than  a  classical 
algorithm  with  the  same  number  of  stages  and  step  length.  The  exact  cost 
of  the  new  integration  scheme  can  be  found  in  [6].  however,  if  we  compare 
the  performance  of  the  geometrically  stable  method  with  the  classical  fourth 
order  Runge-Kutta  algorithm  on  the  problem  of  spherical  elastica,  it  turns 
out  that  a  slightly  increased  step  length  for  the  new  method  suffices  to  out¬ 
perform  the  classical  integration  method.  Using  MATLAB  implementations 
of  both  algorithms  we  found  that  the  geometrically  stable  method  needs  an 
average  value  of  1001  Mflops  for  one  con  plete  step  while  the  fourth  order 
Runge-Kutta  only  uses  an  average  vahn  of  768  Mflops  per  step.  Therefore, 
one  can  statistically  achieve  the  same  performance  by  using  an  increased  step 
length  of  /i  ~  1.3038 /iflK  for  the  new  method.  Taking  into  account  that  the 
proposed  algorithm  not  only  delivers  points  that  lie  exactely  on  the  sphere 
but  that  are  also  approximately  equally  spaced  along  the  tracked  curve,  the 
geometrically  stable  method  seems  to  be  a  good  choice  for  the  integration  of 
spherical  elastica. 

A  performance  comparison  of  the  new  integration  scheme  with  a  more 
sophisticated  classical  method,  the  IMSL  Runge  Kutta  implementation,  can 
be  found  in  [7]. 


5  Classification  of  spherical  elastica 

Acting  on  the  suggestion  of  D.  Bernoulli,  L.  Euler  derived  differential  equa¬ 
tions  for  plane  elastica  and  classified  the  fundamental  forms  of  these  curves 
(see  [9],  [1 1]).  A  curvature  analysis  of  the  various  fundamental  cases  has  been 
given  in  [3]. 

In  this  section  we  classify  the  forms  of  spherical  elastica  based  on  the 
differential  equation  (21)  for  the  geodesic  curvature.  This  equation  is  of  the 
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same  form  as  the  equation  for  the  curvature  of  plane  elastica  and  can  be 
solved  in  terms  of  Jacobi’s  elliptic  functions  in  the  form: 

dnl/xjjd'^  I  ^  ) 

where  the  positive  parameter  P  of  the  elliptic  function  is  given  by 

,3  =  ,29) 

(see  [3]).  The  parameter  represents  the  amplitude  of  the  periodic  curva¬ 
ture  function  and  denotes  the  value  at  which  =  ^m- 

To  obtain  a  representation  of  the  curvature  in  terms  of  Jacobi’s  functions 
with  parameter  P  smaller  than  1 ,  one  uses  the  formula 

^^g(s)  =  Km  Cn(\/(K2^-f  2/r2  -(t)/2{s  -  S,n)  \  (30) 

if  cr  <  0.5k?„  -f  2/r^.  Since  the  function  cn  has  zeros  while  dn  is  positive,  the 
above  case  distinction  reflects  the  main  division  of  elastic  curves  into  those 
where  the  geodesic  curvature  changes  sign  and  the  other  with  constant  sign 
of  their  geodesic  curvature. 

The  change  of  the  forms  of  spherical  elastica  while  Km  is  fixed  and  a  in¬ 
creases  is  shown  by  figures  1-8.  The  maximum  value  of  the  tension  parameter 
cr  for  a  real  elastic  curve  on  a  sphere  is  according  to  (29)  cr  =  k^^  -f-  2/r^. 
This  choice  of  a  corresponds  to  the  dashed  circle  that  is  shown  in  all  figures 
for  the  purpose  of  orientation.  The  second  curve  in  figure  1  ha^  a  negative 
tension  parameter  of  high  absolute  value  (cr  =  —10000).  In  comparison  to 
figure  2  where  a  =  —30  we  observe  that  decreasing  the  tension  parameter 
has  the  effect  of  lowering  the  amplitude  of  the  curve  as  it  is  known  from 
the  euclidean  case  (see  [.3]).  The  oscillation  of  the  curve  in  figure  1  is  of 
too  low  amplitude  to  be  visually  observed  and  the  curve  can  not  be  distin¬ 
guished  from  the  geodesic  determined  by  the  initial  conditions.  (Note,  that 
the  geodesic  curvature  (30)  is  far  from  getting  flat  for  cr  —»  —cx)  but  in  fact 
approaches  a  cos  function  that  oscillates  with  a  period  that  decreases  with 

CT.) 

A  curve  with  positive  tension  parameter  is  shown  in  figure  3  where  cr  = 
2/;’^  =  2.  The  displayed  curve  is  a  special  case  because  the  parameter 
\/P  =  1/2.  Here  the  geodesic  curvature  is  given  by  the  lemniscate  function: 

Kp(5)  =  K„,  coslemn(K„(s  -  s,„)/2). 

In  figures  4  and  5  it  is  illustrated  that  with  increasing  positive  a  the 
bays  of  the  curves  start  to  overlap  until  a  figure-8-configuration  is  reached 
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Figure  1:  =  4,  <7  =  —10000 


where  all  double  points  of  the  curve  coincide.  This  happens  for  a  parameter 
a  «  7.849. 

While  <7  increcises  further  the  curve  proceeds  through  the  figure-8-shape 
forming  a  series  of  loops  with  alternating  sign  of  geodesic  curvature  (see 
figure  6).  These  loops  recede  from  each  other  until  in  the  limiting  case,  when 
<7  =  2/r^  +  0.5k^,  the  curve  forms  a  single  loop  (see  figure  7).  Here  the 
geodesic  curvature  is  given  by 

Kgis)  =  K,nSech(K^{s  -  Sm)/2). 

Figure  8  shows  that  the  single  loop  transforms  into  a  series  of  loops  with 
the  same  sign  of  geodesic  curvature.  With  increasing  o  the  loops  come  closer 
together  and  finally  collapse  into  a  circle  when  a  =  2/r^  +  kJ;,  . 
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